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Summary.  This  research  provides  the  results  of  testing  common  Krylov  subspace  linear  iterative 
solvers  and  preconditioners  in  a  parallel  environment  using  the  Portable,  Extensible  Toolkit  for 
Scientific  Computation  (PETSc)  software  on  extremely  heterogeneous,  unsaturated  flow 
problems.  A  three-dimensional  (3-D)  model  of  a  levee  with  a  root  zone  embedded  at  the  toe 
served  as  the  test  problem,  and  the  runs  were  done  on  the  Cray  XE6  with  128  cores. 


1  INTRODUCTION 

The  record-breaking  floods  of  2011  brought  the  importance  of  the  levee  systems  protecting 
cities  and  farmland  near  our  rivers  to  the  forefront.  Sand  boils  from  underseepage  have  become 
common  terminology  for  many  Americans.  This  research  illustrates  the  importance  of  using  high 
performance,  parallel  computing  to  analyze  levee  performance  during  flood  conditions.  Both 
two-dimensional  (2-D)  and  three-dimensional  (3-D)  numerical  models  of  levees  with  embedded 
tree  root  systems  based  on  the  finite  element  method  were  built  and  run.  Specifically,  an  initial 
analysis  of  the  effect  of  woody  vegetation  on  levees  with  respect  to  the  initiation  of  internal 
erosion  and  the  variability  of  hydraulic  conductivity  was  completed. 
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These  types  of  simulations  create  significant  challenges  when  using  a  time-stepping  algorithm 
such  as  Euler's  implicit  algorithm.  This  is  because  an  ill-conditioned  linear  system  of  equations 
must  be  solved  at  each  nonlinear  iteration  for  each  time-step.  Coincident  to  this  effort,  a  basic 
research  project  was  conducted  to  determine  the  best  Krylov  subspace  iterative  linear  solvers  for 
navigation  locks  and  watershed  simulations1 ,2,3,4.  This  research  adds  to  the  knowledge  base  by 
focusing  on  solvers  for  woody  vegetation  on  levees.  The  Portable,  Extensible  Toolkit  for 
Scientific  Computation  (PETSc)5  library  was  used  in  both  efforts. 

2  LEVEE  MODEL 

Fig.  1  shows  the  cross  section  with  material  types  for  the  Pocket  Levee  in  Sacramento,  CA, 
and  Table  1  gives  the  saturated  hydraulic  conductivities.  Fig.  2  provides  an  example  of  the  2-D 
mesh  that  was  used,  and  Fig.  3  shows  a  3-D  model  of  the  levee  that  was  generated  by  extruding 
the  2-D  cross  section  in  the  third  dimension  multiple  times. 


Levee  sand 


Slurry  wall 


Clay  mixed  with  sand 


Figure  1:  Cross  section  of  the  Pocket  Levee,  Sacramento,  CA 


Material 

kH  (cm/sec) 

kH  (ft/day) 

kv  (cm/sec) 

kv  (ft/day) 

Levee  sand 

8.00  x  10'3 

22.7 

2.00  x  10'3 

5.67 

Clay  and  silty  clay 

8.00  x  10“4 

2.27 

2.00  x  10-4 

0.568 

Clay  mixed  with 
sand 

3.00  x  10'5 

0.085 

1.00  x  10'5 

0.0283 

Aquifer  sand 

8.00  x  10'2 

226.7 

2.00  x  10'2 

56.7 

Gravel 

2.00  x  10'2 

56.7 

2.00  x  10'2 

56.7 

Silt 

1.00  x  10-4 

0.283 

1.00  x  10-4 

0.283 

Slurry  wall 

1.00  x  10'6 

0.00283 

1.00  x  10'6 

0.00283 

Table  1:  Hydraulic  conductivities  used  for  materials  identified  in  the  cross  section  shown  in  Figure  1 
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Figure  2:  Portion  of  finite  element  mesh  of  levee  cross  section 


Figure  3:  Example  of  3-D  mesh  generated  from  2-D  model 
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3  ROOT  ZONE  MODEL  TYPES 
3.1  Changes  in  hydraulic  conductivity 

As  shown  in  Fig.  4,  tree  root  zones  for  the  2-D  analyses  were  modeled  by  creating  near 
rectangular  6-ft  x  5-ft  areas  as  estimated  from  geophysical  surveys6.  Hydraulic  conductivity  for  a 
given  tree  root  zone  was  assigned  to  each  element  using 

^ veg  P^no-veg  (1) 

*  «  4  *  *  « 


Figure  4:  Tree  placement  for  Pocket  Levee,  Sacramento,  CA 

where  kveg  is  the  modified  hydraulic  conductivity,/?  is  a  positive  parameter  with  recommended 
values6  of  0.01  <  J3  <  100 ,  and  kno_ veg  is  the  original  hydraulic  conductivity  without  the  tree. 

3.2  Macropore  heterogeneity 

Another  way  of  modeling  root  systems  is  to  extend  the  root  zone  concept  by  filling  the  root 
zone  with  small  triangular  elements  in  2-D  and  prism  elements  in  3-D  of  approximately  1  inch  in 
each  direction,  but  for  this  model,  the  hydraulic  conductivity  of  each  small  element  is  randomly 
varied.  For  the  ith  small  element  in  the  root  zone,  a  random  number  was  generated,  o<C<i> 

then  /?,  was  computed  for  the  ith  element  using 

Pi  =  104^2  (2) 

Fig.  5  shows  this  concept  as  implemented  in  3-D.  In  fact,  this  dataset  provides  significant 
challenges,  and  it  is  the  one  used  for  testing  to  determine  the  best  linear  solver  and 
preconditioner  for  seepage  in  levees  with  woody  vegetation.  The  mesh  has  3,017,367  nodes. 

4  COMPUTATIONAL  CHALLENGE 

The  computational  challenges  for  this  dataset  are  summarized  as  follows: 
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Figure  5:  Small  3-D  1-inch  elements  in  root  zone 


•  Hydraulic  conductivity  in  the  saturated  flow  region  ( ks )  for  different  soils  (such  as  sand 

and  clay)  can  vary  several  orders  of  magnitude.  The  values  shown  in  Table  1  vary  four 
orders  of  magnitude. 

•  Horizontal  hydraulic  conductivity  is  often  four  times  or  more  greater  than  vertical 
hydraulic  conductivity. 

•  For  flow  in  the  unsaturated  zone,  hydraulic  conductivity  is  modeled  by 

k  =  krks  (2) 

where  kr  is  the  relative  hydraulic  conductivity  and  0  <  kr  <  1 .  kr  varies  several  orders 
of  magnitude. 

•  The  heterogeneous  root  zone  shown  in  Fig.  5  can  have  adjacent  elements  in  which  the 
hydraulic  conductivity  can  differ  by  as  much  as  four  orders  of  magnitude. 

•  All  the  above  points  combine  to  create  ill-conditioned  linear  systems  of  equations  that 
were  solved  at  each  nonlinear  iteration. 
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5  RESULTS 

Four  linear  iterative  solvers  and  four  preconditioners  were  run  using  128  cores  on  the  Cray 
XE6.  The  system  of  equations  is  symmetric  and  positive-definite.  The  solvers  used  are  conjugate 
gradient  (CG),  conjugate  residual  (CR),  bi-CG  stabilized  (BI),  and  Generalized  Minimal 
Residual  (GM).  The  preconditioners  used  are  Jacobi,  block  Jacobi,  Additive  Swartz  Method 
(ASM),  and  Boomer  Algebraic  Multigrid  (AMG).  Fig.  6  shows  a  plot  of  the  timing  results,  and 
Table  2  gives  details  on  iteration  counts  and  running  times.  Runs  were  made  without  (labeled  0) 
and  with  (labeled  1)  the  presence  of  a  tree.  A  dash  indicates  that  the  solution  failed  either  by  not 
converging  after  100,000  iterations  or  the  solver  gave  error  messages. 


□  Jacobi 

□  Block  Jacobi 

□  ASM 

□  Boomer  AMG 


CGO  CGI  CR0  CR1  BIO  BI1  GMO  GM1 


Figure  6:  Running  times  (sec)  for  preconditioner  and  solver  options 


6  CONCLUSIONS 

The  conclusions  drawn  from  the  data  in  Table  2  are  as  follows: 

•  The  presence  of  a  tree  root  made  the  linear  system  of  equations  harder  to  solve  with  some 
solver  and  preconditioner  options  failing  when  a  tree  was  present.  In  particular,  no 
preconditioner  worked  with  CR  when  a  tree  was  present. 

•  The  ASM  preconditioner  did  not  work  with  either  CG  or  CR  for  any  of  the  datasets. 

•  The  Boomer  AMG  preconditioner  performed  significantly  better  than  the  other 
preconditioners. 

•  The  Boomer  AMG  preconditioner  performed  equally  well  with  CG,  BI,  and  GM. 

•  BI  and  GM  took  longer  to  run  than  CG  and  CR. 
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Conjugate  Gradient  (CG) 

Tree  Root 

Preconditioner 

Iterations 

Time  (sec) 

No 

Jacobi 

11,045 

36.5 

Yes 

Jacobi 

- 

- 

No 

Block  Jacobi 

3,204 

17.2 

Yes 

Block  Jacobi 

- 

- 

No 

ASM 

- 

- 

Yes 

ASM 

- 

- 

No 

Boomer  AMG 

51 

5.8 

Yes 

Boomer  AMG 

70 

7.4 

Conjugate  Residual  (CR) 


Bi-CG  Stabilized  (BI) 


No 

Jacobi 

11,152 

73.6 

Yes 

Jacobi 

15,007 

101.0 

No 

Block  Jacobi 

3,034 

33.3 

Yes 

Block  Jacobi 

3,364 

38.3 

No 

ASM 

8,779 

122.2 

Yes 

ASM 

9,365 

130.6 

No 

Boomer  AMG 

28 

6.2 

Yes 

Boomer  AMG 

30 

6.6 

GMRES  (GM) 

No 

Jacobi 

- 

- 

Yes 

Jacobi 

- 

- 

No 

Block  Jacobi 

33,572 

344.0 

Yes 

Block  Jacobi 

33,022 

322.7 

No 

ASM 

16,830 

194.2 

Yes 

ASM 

20,070 

230.2 

No 

Boomer  AMG 

49 

5.9 

Yes 

Boomer  AMG 

51 

6.1 

Table  2:  Performance  of  preconditioners  and  solvers 
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